Evolutionary Model of Species Body Mass Diversification 



o 
o 

(N 

c 

a 

(N 
(N 



cm: 

6 : 
X) : 

(N : 
> ■ 

^t- : 

o ■ 

^ : 

00 ■ 

o : 

oo , 
o ■ 



X 



A. Clauset 1 ^ and S. Redner^'Q 

1 Santa Fe Institute, 1399 Hyde Park Road, Santa Fe NM, 87501, USA 
Center for Polymer Studies and Department of Physics, Boston University, Boston, MA 02215, USA 

We present a quantitative model for the biological evolution of species body masses within large 
groups of related species, e.g., terrestrial mammals, in which body mass M evolves according to 
branching (speciation) , multiplicative diffusion, and an extinction probability that increases logarith- 
mically with mass. We describe this evolution in terms of a convection-diffusion-reaction equation 
for InM. The steady-state behavior is in good agreement with empirical data on recent terres- 
trial mammals, and the time-dependent behavior also agrees with data on extinct mammal species 
between 95 - 50 Myr ago. 
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Animals — both extant and extinct — exhibit an enor- 
mously wide range of body sizes. Among extant ter- 
restrial mammals, the largest is the African savannah 
elephant (Loxodonta africana africana) with a mass of 
10 7 g, while the smallest is Remy's pygmy shrew (Sun- 
cus remyi) at a diminutive 1.8 g. Yet the most probable 
mass is 40 g, roughly the size of the common Pacific rat 
(Rattus exulans), is only a little larger than the small- 
est mass. More generally, empirical surveys suggest that 
such a broad but asymmetric distribution in the number 
of species with adult body mass M typifies many ani- 
mal classes [H, 0, 0, 0, [H , including mammals, birds, fish, 
insects, lizards and possibly dinosaurs. 

What mechanisms cause species mass distributions to 
assume such a shape? A satisfactory answer would have 
wide implications for the evolution and distribution of 
the many other species characteristics that correlate with 
body mass, including life span, metabolic rate, and ex- 
tinction risk 0, 0] . Previous explanations for the species 
mass distribution focused on detailed ecological, environ- 
mental and species-interaction assumptions Q. However, 
empirical data present confusing and often inconsistent 
support for these theories, and none explicitly address 
how species body mass distributions diversify in time. 

In this Letter, we construct a physics-based convection- 
diffusion-reaction model to account for the evolution of 
the species mass distribution. The steady-state behav- 
ior of this model was recently solved to explain the 
species mass distribution for recent terrestrial mammals 
and birds [8, 9], where recent is conventionally defined as 
within the past 50,000 years [lfj. Here we substantially 
extend this approach to give predictions on mammalian 
body mass evolution that are in good agreement with fos- 
sil data. This model can further be used to estimate the 
historical rates of body mass diversification from fossil 
data, which are otherwise estimated using ad hoc tech- 
niques. To illustrate this application, we estimate body 
mass diversification rates from our fossil data, which are 
in good agreement with estimates of genetic diversifica- 
tion from molecular clock methods llj. Although our 
focus is on terrestrial mammal evolution, this model can, 



in principle, be applied to any group of related species. 

Let c(x, t) denote the number (density) of species hav- 
ing logarithmic mass x — InM at a time t; we use x as 
the basic variable in keeping with widespread usage in the 
field 0] . Our model incorporates three fundamental and 
empirically-supported features of biological evolution. 

1. Branching multiplicative diffusion [J 0]: each 
species of mass M produces descendant species 
(cladogenesis) with masses AM, where A is a ran- 
dom variable, and the sign of the average (In A) 
denotes bias toward larger or smaller descendants. 
Empirical evidence [8J1 suggests that (In A) > 
(known as Cope's rule [l|) for terrestrial mammals. 

2. Species become extinct independently, and with a 
probability that increases weakly with mass 12 1 . 



3. No species can be smaller than a mass M m ; n , due, 
for example, to metabolic constraints [l3j |. 

The production of descendant species corresponds to 
growth in the number of species at a rate k that is pro- 
portional to the density c itself. Similarly, the prob- 
ability p{x) that a species of logarithmic mass x be- 
comes extinct may also be represented by a loss term 
that is proportional to c, but with a weak mass depen- 
dence [8(. We make the simple choice of linear depen- 
dence p(x) = A + Bx (but see below). With these in- 
gredients, c(x,t) obeys the convection-diffusion-reaction 
equation in the continuum limit 
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with bias velocity v = (In A) and diffusion coefficient 
D = ((In A) 2 ), and where k — A sets the absolute scale of 
species body mass frequencies. 

To solve Eq. ([I}, we substitute the eigenfunction ex- 
pansion c(x,t) — J2 n A n C n {x) e -7 "*, which yields 
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FIG. 1: (color online) Steady-state solution of the model 
for the species body mass distribution and suitably binned 
empirical data for 4002 recent terrestrial mammals (shown as 
a normalized histogram with 50 logarithmically-spaced bins) . 



where the prime denotes differentiation with respect 
to x, fi = v/D, a = (k - A)/D, and ft = B/D. 
We eliminate the first derivative term by introducing 
ip n = e - ^/ 2 C n and then we use the scaled variable 
z = j3 1 / 3 x - (3~ 2/3 ~ x + D") t0 tran sform Eq. © 
into Airy's differential equation, ip" — zip = for 
each eigenfunction ip n . The general solution is ip(z) = 
ai Ai(z) + a2Bi(z), where Ai(;z) and Bi(z) are the Airy 
functions; here the prime now denotes differentiation 
with respect to z. Since there can be no species with in- 
finite mass and Bi(z) diverges as z — > oo, we set ai = 0, 
while a\ is determined by the initial condition. (We 
could also incorporate a fixed maximum species body 
mass M max , but the analysis is more complicated without 
revealing additional insights.) 

Since the species density vanishes at the minimum 
mass M m j n , the argument of ip n must equal 
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The functions ip n = Ai(z n ) form the complete set of 
states for the eigenfunction expansion [lj| . The first few 
zeros z n are at (roughly) -2.3381, -4.0879, -5.5205 and 
-6.7867 for n = 0, 1, 2,3 14]. (For computing the distri- 
bution, we tabulated the first million zeros numerically 
using standard mathematical software.) The correspond- 
ing decay rates 7„ are then given by j n = D P 2 ^^(zo — z n ), 
with 70 = to give the steady state solution [9j. These 
rates form an increasing sequence so that the higher 
terms in the eigenfunction expansion decay more quickly 
in time. Finally, solving Eq. ([3]) for a and plugging the 



FIG. 2: (color online) Data on species body mass for 569 
North American terrestrial mammals from 95 - 50 Myr 
ago. Each horizontal segment represents a species, and end- 
points denote its first and last appearance in the fossil record. 
The superimposed curve shows the average of In M for these 
species (smoothed with an exponential kernel). 



result into the definition of z, we can eliminate the scale 
parameter a and write z — z n + /3 1 / 3 (a; — x m j n ). Thus 
each eigenfunction has the form 



C n (x) oc e"*/ 2 Ai L + p 1 / 3 ^ - x min ) 



(4) 



in which Ai(z) ~ e _2z ' / 3 for large z. The competition 
between this decay and the prefactor e^/ 2 in C n (x) con- 
tributes to the broadness of the species mass distribution 
and the location of the most probable mass (Fig. |T|) . Par- 
enthetically, the asymptotic decay of the cigenfunctions 
depends only weakly on the form of the extinction prob- 
ability p(x). For instance, if we choose p(x) = A + Bx 5 , 
then as z — > oo, the eigenfunctions decay as e~ z + 1 . 

Suppose that a given group of animals began its evolu- 
tionary history with a single species of mass Mq (initial 
condition c(x,t = 0) = S(x — xq), with xq — InAfo), 
after which speciation occurs according to the dynamics 
of Eq. ((T|). We use the fact that the Airy differential 
equation is a Sturm-Liouville problem [l5[ so that the 
{C n } form a complete and orthonormal set. Following 
the standard prescription to determine the coefficients 
of the eigenfunction expansion, the full time-dependent 
solution is 



c(x, t) — 



(x )C n (x) 



(5) 



where each C n {x) is given by Eq. In the long-time 
limit, all the decaying eigenmodes with n > become 
negligible and the species mass distribution reduces to 
c(x,t — > oo) oc Cq (x) . 
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FIG. 3: (color online) A comparison of our model predictions 
(dashed) from Eq. ((5} and the mass distributions of extinct 
North American mammal species (solid) in nine consecutive 
time ranges. (To get reasonable results with sparse empirical 
data, we smooth the distributions with a Gaussian kernel.) 

We test our predictions for the species body mass dis- 
tribution by comparing with available empirical data. In 
the steady state, our model is characterized by three 
parameters: M m i n , the mass of the smallest animal, /i, 
which controls the tendency of descendant species to be 
larger or smaller than their ancestors, and 0, which con- 
trols the dependence of the extinction rate on mass. The 
former two can be estimated directly from fossil data, 
while the latter is typically estimated by matching the 
steady-state solution to the recent data. For terrestrial 
mammals, we previously found /i w 0.2, (3 rs 0.08, while 
Afmin ~ 2 g [9] . Using these parameter values in the long- 
time limit, we obtain a good agreement between the pre- 
dictions of the model and the species mass distribution 
of recent terrestrial mammals [10| (Fig. Q]) . 

Our model also makes predictions about the way the 
body mass distribution changes over time, which can be 
tested with fossil data. Drawing on data from the best 
available source on the evolution of mammalian body 
masses [l6[ , we plot in Fig. [2] the durations and body 
masses of 569 extinct mammal species from 95 - 50 Myr 
ago. This period includes the Cretaceous-Paleogene (KP) 
boundary of 65.5 Myr ago that marks a mass extinc- 
tion event during which more than 50% of then extant 
species became extinct, including non-avian dinosaurs — 
the dominant fauna for the preceding 160 Myr — and is 
the subject of many studies regarding the diversification 
of mammals. (The Cretaceous period is conventionally 
abbreviated "K," after the German translation Kreide.) 

Using the same model parameters as above, and setting 
Mq = 2g, the estimated size of the first mammal fl7j j . 
Fig. [3] shows good agreement between model predictions 
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FIG. 4: Estimated relation between model time to histori- 
cal time, indicating that the broadening of the mammalian 
species mass distribution began approximately 10 Myr before 
the KP boundary when non-avian dinosaurs became extinct. 

from Eq. (O and empirical data from Fig. [5] To sim- 
plify the comparison, we divided the period from 95 - 50 
Myr ago into nine bins of equal durations, and tabulated 
the distribution of species extant during each of these 
bins. In each of the nine panels of Fig. [3l we give both 
the historical time period and the corresponding model 
time that yields a good qualitative match to the data. 
(The fitting can be made more objective using standard 
techniques, but the results are largely the same.) 

The relation between historical time and model time is 
itself interesting (Fig. H]) . During the first 20 Myr (95 - 
75 Myr ago), the species mass distribution is almost sta- 
tionary, with model time advancing by only At = 2.30; 
curiously between 90 - 80 Myr ago, model time does 
not advance at all. This period of near stasis may indi- 
cate a lull in the evolutionary dynamics, perhaps due to 
implicit competition from larger species, e.g., dinosaurs. 
Over the next 20 Myr (75 - 55 Myr ago), however, the 
distribution broadens considerably (model time advanc- 
ing by At = 10.74), and comes to closely resemble the 
recent distribution (Fig. [I]). This correspondence sug- 
gests that the diversification of mammalian body masses 
into their recent state began at least 75 Myr ago, roughly 
10 Myr before the KP boundary and the extinction of 
the non-avian dinosaurs. This estimate of the timing 
of body mass diversification for mammals agrees closely 
with some estimates of the timing of mammalian ge- 
netic diversification from studies of molecular clocks [ijj , 
and supports the notion that mammals were diversifying 
prior to and independently of the KP boundary itself [l8| . 
Whether these two forms of diversification are causally 
linked, however, is unknown. 

Here we have held all model parameters constant while 
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adjusting model time to fit the data. In principle, how- 
ever, model time could advance steadily while varying 
some model parameters, perhaps to reflect large-scale 
changes or trends in the selection pressures on species 
body size. Empirical evidence supports a stable value of 
M m i n (Fig. [2]), but little is currently known about how 
or why fi and (3 may have varied. 

As is typically the case with historical inference using 
the fossil record, a few caveats are in order. Our fossil 
data are derived from the well-studied North American 
region using modern dental techniques, which are less 
prone to biases than older techniques. Still, some biases 
and sampling gaps likely persist, and these may explain 
the slight overabundance of large species, and under- 
abundance of small species, in more recent times (Fig. [3]). 
More significantly, recent fossil discoveries suggest that, 
since mammals originated 195 Myr ago 17| . mammalian 
diversification has proceeded in several waves, and the 
vast majority of species groups from the earlier waves 
are now extinct. Our data cover only the most recent di- 
versification, in which therian (placental and marsupial) 
mammals largely replaced the then dominant non-therian 
mammal groups [191 ]. Unfortunately, suitable data on 
these waves of diversification is not currently available, 
and the data we do have is sparse in its coverage of non- 
therians. Thus, the application of our model to infer 
diversification rates should be considered as a proof-of- 
concept, illustrating that a physics-style model can shed 
considerable light on evolutionary dynamics by placing 
the fossil record within a theoretical framework. 

In summary, the broad distribution of body masses 
for mammals appears to be well described by a simple 
convection-diffusion-reaction model that incorporates a 
small number of evolutionary features and constraints. 
Indeed, our model does not account for many canoni- 
cal ecological and microevolutionary factors, such as in- 
terspecific competition, geography, predation, and pop- 
ulation dynamics [3j. The fact that our model agrees 
with species mass data suggests that the contributions of 
the above-mentioned processes to the global character of 
body mass distributions can be compactly summarized 
by the parameters /i and (3 in our model. Despite the 
crudeness of the model, the agreement between its pre- 
dictions and the available empirical data is satisfying. 

Our model opens up intriguing directions for theo- 
retical descriptions of evolutionary dynamics. For in- 
stance, our model ignores the populations of individuals 
within each species; estimating the sizes of these pop- 
ulations from body mass vs. population density scaling 
relationships |20l| and body mass vs. home-range size re- 
lationships [21| may allow us to calculate both the to- 
tal biomass contained in a group of related species, and 
its temporal dynamics during diversification. Similarly, 
when paired with scaling relations between body mass 
and metabolism [l3j], we may be able to calculate the 



total metabolic flux of a taxonomic group. 
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